#!/usr/bin/R

#function that determines flow rates across 'boxs' (weirs/pumps)
###############################################################################################
find.flow = function(my.h, my.box) { 
	my.type=as.integer(box[,1])
	my.weir=as.numeric(box[,2])
	my.width=as.numeric(box[,4])
	my.roughness=as.numeric(box[,5])
	my.pump=as.numeric(box[,6])
	my.upstream=as.integer(box[,7])
	my.downstream=as.integer(box[,8])
		
        #initialize local variables
        my.q=array(0,dim=n.boxs)
        my.inflow=array(0,dim=length(my.h))
        my.outflow=array(0,dim=length(my.h))

        for (b in 1:n.boxs) {
		
		upstream=my.upstream[b]
		downstream=my.downstream[b]
 
		if (my.type[b]==2) { #box is a pump
                        my.q[b]=my.pump[b]
                }
                
		if (my.type[b]==1) { #box is a weir
                        if (my.h[upstream]>my.weir[b] & my.h[downstream]<=my.weir[b]) {
                                #forward flow, not submerged
                                my.q[b]=3.33*my.width[b]*(my.h[upstream]-my.weir[b])^(3/2)
				#print(paste(b,'forward',my.q[b]))
                        }
                        if (my.h[upstream]<=my.weir[b] & my.h[downstream]>my.weir[b]) {
                                #backward flow, not submerged
                                my.q[b]=(-1)*3.33*my.width[b]*(my.h[downstream]-my.weir[b])^(3/2)
				#print(paste(b,'backward',my.q[b]))
                        }
                        if (my.h[upstream]>my.weir[b] & my.h[downstream]>my.weir[b]){
                                #submerged sharp-crested weir formula
                                my.q[b]=abs(3.33*my.width[b]*sqrt(abs(my.h[upstream]-my.h[downstream]))*
                                        ((my.h[upstream]-my.weir[b])+(my.h[downstream]-my.weir[b])*.381))
                                #adjust sign for direction
                                if (my.h[upstream]<my.h[downstream]) {
                                        my.q[b]=my.q[b]*(-1)
                                }
				#print(paste(b,'submerged',my.q[b]))
                        }
                }

                my.q[b] = my.q[b]*my.roughness[b] #multiply flow by empirical roughness factor
		my.outflow[upstream]=my.outflow[upstream]+my.q[b]
                my.inflow[downstream]=my.inflow[downstream]+my.q[b]
        }
        return(list(my.inflow,my.outflow))
} #end function find.flow()
###############################################################################################


